Проект №4

Портфельная оптимизация

Маргарита Бурова, Артем Дремов, Камель Есимова, Юрий Кульчицкий, Ольга Сабирова

ММОС16

Задача 1. Quadratic programming

Постановка задачи

Решить задачу квадратичного программирования $$ \begin{align}{c} \min -\bar{p}^Tx + \mu x^T \Sigma x\\ \sum x_i = 1,\ x_i \geq 0 \end{align} $$ для большого количества различных положительных $\mu$ (например, 100 значений, логарифмически распределенных от 1 до $10^7$). Нарисовать графики $\bar{p}^T x$ от стандартного отклонения $(x^T \sigma x)^{1/2}$, а также картинки, как на (рис. 4.12, стр. 187, Boyd).

Первый набор данных

$$\overline{p}=\begin{bmatrix} 0.12\\ 0.10\\ 0.07\\ 0.03 \end{bmatrix}, \ \ \ \sum= \begin{bmatrix} 0.0064 & 0.0008 & -0.0011 & 0\\ 0.0008 & 0.0025 & 0 & 0\\ -0.0011 & 0 & 0.0004 & 0\\ 0 & 0 & 0 & 0 \end{bmatrix}$$

В наборе четыре инструмента, один из которых безрисковый и дает доход в 3%.

Методы решения задач квадратичного программирования

Барьерный метод: описание

Для использования этого метода задача переформулируется следующим образом: $$ \begin{array}{rl} \text{minimize} & (m / \varepsilon) f_0(x) + \phi(x)\\ \text{subject to} & Ax = b \end{array} $$

Идея метода состоит в решении последовательности минимизационных подзадач: вычисляется значение $x^*(t)$ для возрастающей последовательности из $t$, пока $t < m / \varepsilon$.

Барьерный метод: алгоритм

Алгоритм выглядит следующим образом: сначала задаются $x, t := t^{(0)} = 0, \mu > 1, \varepsilon > 0$. Затем повторяются следующие операции:

1. Центрируется шаг. Вычисляется $x^*(t)$, минимизируя $tf_0 + \phi$ относительно $Ax = b$, начиная с $x$.
2. Обновляется значение $x := x^*(t)$.
3. Проверяется выполнение критерия остановки: выход из цикла при $m/t < \varepsilon$.
4. Увеличивается $t$: $t := \mu t$.

• Метод останавливается, когда $f_0(x) - p^* \leq \varepsilon$ (следует из $f_0(x^(t)) - p^ \leq m / t))

• Центрирование обычно делается методом Ньютона, начиная от текущего значения $x$.

• Выбор $\mu$ подразумевает компромисс: большое значение $\mu$ означает меньшее число внешних итераций, большее -- большее; обычно выбирают значения от 10 до 20.

• Используется несколько эвристик для выбора $t^{(0)}$.

Барьерный метод: анализ сходимости

Число внешних итераций (центрирования) задается, как $$ \bigg\lceil \dfrac{\log(m / \varepsilon t^{(0)}))}{\log \mu}\bigg\rceil $$ плюс изначальный шаг центрирования (для вычисления $x^*(t^{(0)})$).

Задача центрирования: $$ \min t f_0(x) + \phi(x) $$

Для этого необходимо иметь результаты анализа сходимости метода Ньютона.

• $tf_0 + \phi$ должен иметь замкнутые подуровневые множества для $t \geq t^{(0)}$
• классический анализ предполагает строгую выпуклость и выполнение условия Липшица
• анализ через самосогласованность предполагает самосогласованность функции

Метод Ньютона

Итеративный метод поиска корней уравнения $f'(x) = 0$, при этом $f$ должна быть дифференцируемой. $$ x_{n+1} = x_n - \dfrac{f'(x_n)}{f''(x_n)} $$

Метод внутренней точки

1. Выбирается произвольная точка $x_0$ внутри многогранника решений
2. Находится преобразование $\phi$, для которого $c^T \phi(x) < c^T x$
3. На каждой итерации переходим от $x^k$ к $x^{k+1} = \phi(x^k)$, оставаясь внутри многогранника
4. Когда значение $c^T \phi(x^k)$ станет достаточно маленьким, перейдем из текущей точки в ближайшую вершину многогранника, которая и будет решенем

Результаты численных экспериментов на малом наборе данных

Результаты численных экспериментов на большом наборе данных

Задача 2. Semidefinite programming

Постановка задачи:

Предполагается, что дана нормально распределенная случайная величина $p\sim\mathcal{N}(\bar{p}, \Sigma)$. Требуется сформулировать задачу $$ \begin{array}{c} \max \bar{p}^T x\\ prob(p^Tx\leq 0) \leq \eta\\ \sum x_i = 1,\ x_i \geq 0 \end{array} $$ как задачу выпуклой оптимизации с параметром $\eta < 1/2$, и решить ее для большого числа разных $\eta$ от $10^{-4}$ до $10^{-1}$, нарисовав оптимальные $\bar{p}^T x$ от $\eta$.

$$ \begin{array}{c} \max \bar{p}^T x\\ prob(p^Tx\leq 0) \leq \eta\\ \sum x_i = 1,\ x_i \geq 0 \end{array} $$

Но $$ prob(p^Tx\leq 0) = \Phi\left( \dfrac{-\bar{p}^T x}{\lVert \Sigma^{1/2} x \rVert_2} \right) = 1 - \Phi\left( \dfrac{\bar{p}^T x}{\lVert \Sigma^{1/2} x \rVert_2} \right) $$ Теперь $$ \begin{array}{c} \max \bar{p}^T x\\ -\bar{p}^Tx + \Phi^{-1}(1 - \eta)\lVert \Sigma^{1/2} x \rVert_2 \leq 0\\ \sum x_i = 1,\ x_i \geq 0 \end{array} $$

Second-order cone program (SOCP)

То, что мы вывели раньше, формулируется в терминах SOCP: $$ \begin{array}{rl} \text{minimize} & f^T x\\ \text{subject to} & -\lVert A_i x + b_i \rVert_2 \leq c_i^T x + d_i, i = 1,\ldots,m\\ \end{array} $$ Сформулируем нашу задачу в этих терминах $$ \begin{array}{rl} \text{minimize} & -\bar{p}^T x\\ \text{subject to} & \lVert \Sigma^{1/2} x \rVert_2 \leq \dfrac{1}{\Phi^{-1}(1 - \eta)}\bar{p}^Tx \\ & \sum x_i = 1\\ & \mathbf{diag}(x_i) \succeq 0 \end{array} $$ Ограничения встроились бы в виде стандартной комбинации двух неравенств, которую мы бы преобразовали так, чтобы все выражалось через норму, но здесь этого делать не будем, так как мы найдем новую эквивалентную задачу.

Semidefinite programming (SDP)

$$ \begin{array}{rl} \text{minimize} & c^T x\\ \text{subject to} & x_1 F_1 + x_2 F_2 + \ldots + x_n F_n + G \preceq 0\\ & Ax = b \end{array} $$

Задача такого вида, где $F_i, G \in \mathbb{S}^k$, называется задачей полуопределенного программирования, и именно ее мы и будем решать численно. Для задачи конического программирования второго порядка $$ \begin{array}{rl} \text{minimize} & f^T x\\ \text{subject to} & \lVert A_i x + b_i \rVert_2 \leq c_i^T x + d_i, i = 1,\ldots,m\\ \end{array} $$ существует эквивалентная задача полуопределенного программирования: $$ \begin{array}{rl} \text{minimize} & f^T x\\ \text{subject to} & \left[ \begin{array}{cc} (c_i^T x + d_i)I & A_ix+b_i\\ (A_i x + b_i)^T & c_i^T x + d_i \end{array} \right] \succeq 0,\ i = 1,\ldots,m \end{array} $$ В нее уже можно будет внести наши ограничения явно.

Приведение нашей задачи к виду SDP

$$ \begin{array}{rl} \text{minimize} & -\bar{p}^T x\\ \text{subject to} & \left[ \begin{array}{ccc} \left(\dfrac{1}{\Phi^{-1}(1 - \eta)}\bar{p}^Tx\right)I & \Sigma^{1/2}& 0\\ {\Sigma^{1/2}}^T & \dfrac{1}{\Phi^{-1}(1 - \eta)}\bar{p}^Tx & 0\\ 0 & 0 & 1^T(x) \end{array} \right] \succeq 0,\ i = 1,\ldots,m \end{array} $$

Ограничение на единичную сумму иксов введем, разложив покомпонентно по иксам приведенную выше матрицу, а затем приписав к каждой из матриц справа снизу число 1, и определив соответствующий элемент матрицы $G$, как 1.

Как мы будем разделять эту матрицу на компоненты? Подставляя в нее в качестве $x$ элементы базиса (векторы с единственной ненулевой координатой, равной единице).

Численное решение задачи

Основные методы численного решения задач полуопределенного программирования

• Методы внутренней точки (CSDP, SeDuMi, SDPT3, DSDP, SDPA)
• Методы первого порядка (коническая оптимизация), которые, по идее, более правильные за счет отсутствия необходимости хранить гессиан, но менее точные.
• ConicBundle — пакет, который формулирует задачу в виде задачи негладкой оптимизации и решает ее методом спектральной свертки.
• Алгоритмы, решающие задачу SDP, предполагая, что она является задачей нелинейного программирования.

Мы использовали метод CSDP.

library(Rcsdp)
library(expm)
library(plotly)
p_bar <- c(0.12, 0.10, 0.07, 0.03)
Sigma <- matrix(
  c(
    0.0064, 0.0008, -0.0011, 0,
    0.0008, 0.0025, 0, 0,
    -0.0011, 0, 0.0004, 0,
    0, 0, 0, 0
  ),
  nrow = 4
)

Sigma_sqrt <- sqrtm(Sigma)

optimal_res <- sapply(
  seq(0.0001, 0.1, 0.0001),↔
  function(eta) {
        p_bar.mat <- matrix(p_bar)
    
    C <- list(
      matrix(0, length(p_bar) + 1, length(p_bar)+ 1),
      c(rep(0, length(p_bar)), 1, -1)
    )
    A <- lapply(↔
    )
    b <- (-1) * p_bar
    K <- list(type = c("s", "l"), size = c(length(p_bar) + 1, length(p_bar) + 2))
    sol <- csdp(C, A, b, K, control = csdp.control())
    
    t(p_bar.mat) %*% sol$y
    
  }
)

Error in library(plotly): there is no package called ‘plotly’
Traceback:

1. library(plotly)
2. stop(txt, domain = NA)

df <- data.frame(id = 1:length(optimal_res), return = optimal_res)
plot_ly(df, x = ~id, y = ~return, type="scatter", mode = "lines", height=500)

real.rets <- read.csv(
  "~/Dropbox/Documents/Studies/MA (NRU HSE)/Year 1/Optimization Methods/Project/reveal_Opt/futures_returns.csv",
  sep = ";", check.names = F, stringsAsFactors = F
)

real.rets <- real.rets[real.rets$Date > "2005-01-01", ]
real.rets <- real.rets[, -1]
real.rets <- real.rets
p_bar <- colMeans(real.rets)
real.rets <- real.rets[, which(is.finite(p_bar))]
p_bar <- colMeans(real.rets)
Sigma <- cov(real.rets)
Sigma_sqrt <- sqrtm(Sigma)
p_bar <- p_bar * 100
Sigma_sqrt <- Sigma_sqrt * 100


optimal_res <- sapply(
  seq(0.0001, 0.1, 0.001),
  function(eta) {
        p_bar.mat <- matrix(p_bar)
    
    C <- list(
      matrix(0, length(p_bar) + 1, length(p_bar)+ 1),
      c(rep(0, length(p_bar)), 1, -1)
    )
    A <- lapply(
      1:length(p_bar),
      function(i) {
        x <- rep(0, length(p_bar))
        x[i] <- 1
        x <- matrix(x)
        
        
        this.mat <- cbind(
          rbind(
            as.numeric(1 / qnorm(1 - eta) * t(p_bar.mat) %*% x) * diag(length(p_bar)),
            t(Sigma_sqrt %*% x)
          ),
          rbind(
            Sigma_sqrt %*% x,
            1 / qnorm(1 - eta) * t(p_bar.mat) %*% x
          )
        )
        
        this.vec <- rep(0, length(p_bar) + 2)
        this.vec[i] <- 1
        this.vec[length(p_bar) + 1:2] <- c(1, -1)
        list(
          this.mat,
          this.vec
        )
      }
    )
    b <- (-1) * p_bar
    K <- list(type = c("s", "l"), size = c(length(p_bar) + 1, length(p_bar) + 2))
    sol <- csdp(C, A, b, K, control = csdp.control())
    
    t(p_bar.mat) %*% sol$y
    
  }
)

df <- data.frame(id = 1:length(optimal_res), return = optimal_res)
plot_ly(df, x = ~id, y = ~return, type="scatter", mode = "lines", height=500)

Спасибо за внимание!